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Abstract. An astrophysical population of supermassive black hole binaries is 

thought to be the strongest source of gravitational waves in the frequency range 
covered by Pulsar Timing Arrays (PTAs). A potential cause for concern is that the 
standard cross-correlation method used in PTA data analysis assumes that the signals 
are isotropically distributed and Gaussian random, while the signals from a black hole 
population are likely to be anisotropic and deterministic. Here we argue that while the 
conventional analysis is not optimal, it is not hopeless either, as the standard Hellings- 
Downs correlation curve turns out to hold for point sources, and the small effective 
number of signal samples blurs the distinction between Gaussian and deterministic 
signals. Possible improvements to the standard cross-correlation analysis that account 
for the anisotropy of the signal are discussed. 
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1. Introduction 

The most promising source of signals in the frequency range covered by Pulsar Timing 
Arrays (PTAs) is from a population of supermassive black hole binaries, dominated by 
systems with masses in the range 3 x 10 ''Mq — 3 x IO^Mq, and times to merger in the 
range 10^ years — )■ 10^ years. It had been assumed that the number density of sources as 
a function of frequency, dN/df, would be sufficiently large that the central limit theorem 
would come into play, and that the combined signal would be Gaussian distributed and 
isotropic. However, recent studies based on more realistic population synthesis models 
have shown that the signal is likely to be dominated by a small number of relatively 
nearby sources [H El |3l H] , and as a result, will be non-Gaussian and anisotropic |5l |6]. 
This is a concern since the standard analysis techniques [H [HI E] are based on the 
assumption that the signal is isotropic and Gaussian. 

Here we show that, while the standard approach may not be optimal, it is able 
to detect the signals from isolated black holes, and by extension, populations of black 
holes no matter how sparse. What makes this possible is the rather surprising result that 
the Hellings-Downs correlation curve [7j, which was originally derived for un-polarized, 
isotropic backgrounds, continues to be valid for polarized point sources! Like many 
results that are surprising initially, after a little thought this result starts to make sense 
(it is basically a reflection of the quadrupole nature of the signal), and very soon the 
result becomes obvious, and a soon after that, something everyone knew already. 

While the standard cross-correlation analysis technique can be used to detect the 
signals from a sparse black hole background, it will not be optimal. We consider a 
variety of alternative analysis techniques that may be more effective, and suggest a new 
cross-correlation technique that accounts for the anisotropy of the signal. 

2. Detector Response 

The timing residuals for a Pulsar located in the h — )■ {9p, (pp) direction, induced by a 
plane gravitational wave from a source in the {9, (p) direction, can be expressed as 



where ifj is the polarization angle of the wave relative to the frame deflned by the basis 
vectors u, v that span the plane perpendicular to the propagation direction k, where 



r = - 







u = cos 6 cos (j)x + cos 6 sin (f)y — sin 6 z , 
V = sin (f)x — cos (p y ■ 
The antenna bean pattern functions have the form 



(2) 




F 



X 



1 + k ■ h 
2{u ■ n){v ■ h) 



(3) 



1 + k ■ fi 



Black Hole Backgrounds 



3 



The terms -R+,x are expressed in terms of the anti-derivative, i?+,x of the gravitational 
wave strain h^^^: 

R+^^=H+^^{t)-H+^^{t-L{l + k-h)), (4) 

where L is the distance to the Pulsar from Earth. The two terms in the above equation 
are referred to as the "Earth term" and the "Pulsar term", respectively. For nearby 
sources the plane wave approximation may need to be augmented by the leading order 
spherical wavefront corrections of order L/D, where D is the distance to the source: 

i?+,x = i^+,x (t) - i^+,x (t - L(l + A; ■ n) + ^(1 - (A; ■ hf)) . (5) 
The antenna patterns can be re-written in the alternative, simpler form 
F+ = (l + cos/3) cos(2a) 

= (l + cos/3)sin(2a), (6) 

where (3 = arccos(— A; ■ n) is the angle between the source and the Pulsar, and 
a = arctan(('0 ■ n)/{u ■ h)) is the angle the Pulsar direction makes relative to the u,v 
polarization frame. The timing residuals then take the form 

r = ^ cos(2V^ + 2a) + i?x sm{2ilj + 2a)) (1 + cos /3) . (7) 
3. Correlation Analysis 

The cross-correlation of the timing residuals from two Pulsars can be written as 

{riVj) = ^{{RD cos(2^ + 2ai) cos(2^/' + 2a,) 

+ (Rl) sm{24j + 2ai) sin(2^/' + 2aj)){l + cos/3i)(l + cos/3j), (8) 
where the angle brackets denote the inner product 

{hh2) = JdhJ dt2 hi{t^)K{ti,t2)h2{t2) . (9) 

For stationary signals, the convolution kernel is a function of the lag \ti — ^2!; and the 
inner product can be re-wrriten in the Fourier domain in the familiar form 

(/,.M^f ^"'■<^'''«y/<^'''-'^>' #. (10) 

In ([s]) it has been assumed that (R^R^) = 0, which holds for cosmological stochastic 
backgrounds and binary systems. For isotropic gravitational wave backgrounds it makes 
sense to average the cross correlation over the sky: 

^ I {ur,) dn = {H') a,, , (11) 

where {H"^) = {R^ + R^), and the Hellings- Downs correlation curve aij = a{6ij) is given 
as a function of of the angle 9ij = fi between the Pulsars: 

, , 1-cos/i, /l-cos/i\ 1-cos/i 1,^ ^, ,, 

= ^ In [ ^) ^ + 3 + ^^^^^ ' ^^2^ 
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Figure 1. The auto-correlated signal power {r'^{9p, (pp)) for a single realization of a 
black hole binary population model. The left panel has the "Pulsar terms" turned off, 
while the right panel shows the full response. The Pulsar term adds noise, effectively 
multiplying the Earth-term sky map by 2(1 — cos{6)), where 5 is a random phase. 
Generating the full response at higher angular resolution and then applying a Gaussian 
smoothing yields a sky map nearly identical to the map with the Pulsar term turned off. 
Either way, the anisotropy of the signal is clear, with Pulsars in certain sky directions 
receiving significantly larger signal power than others. 



The delta function - defined such that 6{0) = 1, and is otherwise zero - comes from the 
Pulsar term, which averages to a non-zero value in the auto-correlation. 

For anisotropic signals, such as those produced by a single black hole binary, sky 
averaging is not justified, and the correlation {riVj) will depend on the sky location 
of the source {0,(f)), and the orbital orientation given in terms of the inclination and 
polarization angles {t,4')- It had been assumed that an astrophysical population of 
binaries would combine to yield an isotropic, stochastic background, but this turns out 
not be be the case. Instead the combined signal is dominated by a handful of nearby, 
bright sources, and as shown in Figure 1, the resulting background is highly anisotropic. 
The BH population model used to generate Figure 1 was derived by extracting catalogues 
of merging massive galaxies from the Bertone et al. [10] semianalytic model built on top 
of the Millennium Run [TT]. Galaxies were then populated with super massive black 
holes correlating with the bulge velocity dispersion as given by Gultekin et al. [12]. 
The black holes accrete gas prior to final coalescence and all binaries are assumed to be 
circular and driven by GW emission only in the frequency band relevant to PTA. All the 
steps of the procedure followed to construct the population are given in Sesana et al. [2] . 
The anisotropy seen in Figure 1 is even more pronounced if the signal is broken out by 
frequency bins, where a single source often dominates in a particular bin. The question 
then becomes, what is the best technique to detect such a signal, given that it is neither 
isotropic nor Gaussian? These assumptions underpin the standard analysis techniques 
in both the Frequentist and Bayesian implementations. The Frequentist approach is 
based on the matched filter detection statistic [S] 

p = EE(^^^.)«(%)- (13) 

This statistic is often shifted to have zero mean and scaled to have variance l/iVpairs, 
where A^pairs = A^(A^ — 1) are the number of pulsar pairs. The key idea is that 
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the pairwise correlations are summed together after being multiphed by the expected 
correlation function, which acts like a matched filter. In the Bayesian approach, the 
correlation function enters into the definition of the multi-variate Gaussian likelihood 
function [131 M > 

p = Aexp(^-J2l inf* + f:f,)Cr^df -^J ln(detC)d/j , (14) 

where A is an overall normalization constant and 

CM = SH{f)a.j + SnXm,. (15) 

Here Snif) is the power spectral density of the signal and SnXf) is the power spectral 
density of the noise in the i^^ Pulsar. In the weak-signal limit, Sm ^ Sh, the likelihood 



(14) can be approximated as p = y4'exp(p/2), drawing out the close connection between 



the two approaches. 

4. Isolated Black Hole Binaries 

Before discussing alternative analysis techniques that may be better suited to detecting 
anisotropic, non-Gaussian signals, it is interesting to consider how the standard analysis 
might perform at detecting the signal from an isolated Black Hole binary. To set the 
stage, let us consider the correlations produced in a Pulsar Timing Array with 100 
randomly distributed Pulsars by (i) a single black hole binary; and (ii) an isotropic 
background. To make the comparison equitable, the isotropic signal was restricted to a 
single frequency bin. In Figure 2 the correlations are shown both with and without the 
Pulsar term, and un-binned and binned in the angular separation between the Pulsars. 
The signal strength in each case was scaled to give unit correlation at zero degree 
separation. The results in both cases are very similar. The un-binned correlations show 
significant scatter, while the binned correlations follow the Hellings-Downs correlation 
curve. At first sight it may seem strange that an isolated black hole binary produces a 
correlation pattern that is identical to that produced by an isotropic background, but on 
reflection the result is not surprising. The Hellings-Downs curve is simply a consequence 
of the quadrupolar nature of gravitational waves. In binning the correlations as a 



function of the Pulsar angular separation we are replacing the sky average (11) by 
an average over the Pulsar locations, which in the limit of a large number of Pulsar 
pairs goes over to the integral 

^(/^) ~ /A f {''^i'^j) (^(cos fi — fii ■ hj) dflidQj . (16) 



(47r)2 . 

The Dirac-delta function can be taken care of by adopting a coordinate system where 
the j-Pulsar has coordinates 

Xj = cos (pi (cos 9i sin /i cos A + sin 9i cos /i) — sin (pi sin /i sin A 
Hj = sin 0i(cos 6i sin fi cos A + sin 9i cos /i) + cos (pi sin /i sin A 
Zj = cos 9i cos /i — sin 6i sin /i cos A , (17) 
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which ensure that hi ■ hj = cos yU. Completing the integration over X,(j)i,6i yields 

7(/i) = (i/')«(^). (18) 

Thus, a single black hole produces an identical angular correlation pattern as an isotropic 
stochastic background. Note that the final result is independent of the black hole 
orientation or sky position. Again, this is not surprising since we have integrated the 
Pulsar locations over the celestial sphere, which is equivalent to actively rotating the 
black hole reference frame while holding the Pulsars locations fixed. 
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Figure 2. Simulated noise free correlation curves for (i) an isolated BH binary (ii) 
an isotropic background restricted to a single frequency bin. The top panels show the 
correlations as a function of the angle between the Pulsar pairs without the "Pulsar 
term" ; the middle panels show the full correlation; while the lower panels show the full 
correlations averaged into 5 degree bins. 



5. Astrophysical Black Hole Populations 



In Ref. ^ the applicability of the standard analysis techniques based on (13) and (14) 
for detecting the signals from an astrophysical population of black holes was discussed. 
There the focus was on the non-Gaussianity of the signal, rather than the anisotropy. It 
was noted that the correlations between Pulsars followed the Hellings-Downs correlation 
curve upon averaging over ~ 100 realizations (unsurprising given that the averaging 
restores isotropy), but this result has little practical relevance given that we only get 
to see a single realization. On the other hand, the fact that a single black hole binary 
yields the Hellings-Downs curve means that the standard analysis techniques will be 
effective (though not necessarily optimal) at detecting the signal from a population of 
black holes. And while it is possible to theoretically establish the non-Gaussianity of 
the signal using hundreds of realizations of the population catalogs, it will be difficult 
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to establish in practice with the handful of frequency bins available for the analysis. 
Indeed, the departure from Gaussianity will likely be established by the detection of 
one or more of the brightest signals using single source analysis techniques [Ml [15] . The 
importance of their being few effective samples in the data is illustrated in Figure 3, 
where correlation curves for various simulated signals are shown based on a ten year 
observation period (noise was not added to the signals in Figure 3 so as not to obscure the 
intrinsic scatter from the Pulsar term, but a noise spectrum was used when computing 
the inner products). In these simulations the Pulsar timing noise was assumed to have 
a white spectrum above 6 nHz, and a red spectrum at lower frequencies [16]: 

Snif) = const, (l + (//6nHz)-2) . (19) 
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Figure 3. Simulated noise free correlation curves for (i) an isotropic, White Gaussian 
background using 100 frequency bins (ii) an isotropic Red Gaussian background with 
the spectrum predicted for a BH population (iii) a BH population model. The upper 
panels are raw scatter plots as a function of the angle between the Pulsars, while the 
lower panels average the correlations into 5 degree bins. 



The first panel in Figure 3 shows the correlation curve for an isotropic stochastic 
background with a white spectrum that covers 100 frequency bins. The second panel 
shows the correlation curve for an isotropic stochastic background with a red spectrum 
where the slope was chosen to match that from a population of black hole binaries 
(Snif) ~ f^^^^^)- The third panel shows the correlation curve for a realization of the 
BH population model used to generate Figure 1. Remarkably, the scatter from the BH 
population is less than for an Gaussian stochastic background, as can be seen in the 
histograms of {riVj)/ {H"^) — <y{Oij) shown in Figure 4. 

Having established that the standard correlation analysis is capable of detecting the 
anisotropic, deterministic signal from an astrophysical population of black hole binaries, 
it is worth considering how the analysis can be improved. What we are seeking is an 
analysis technique that has an optimal balance between fidelity in the signal model and 
parsimony in terms of dimensionality. High dimensional models can achieve high fidelity, 
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Figure 4. Histograms of the scatter in the correlation about the HeUings-Downs model 
for an isotropic, red Gaussian background (red, solid line), and for an astrophysical 
population of black hole binaries (blue, dashed line) . 

but at the cost of a larger trials factor (in a Frequentist setting) or Occam factor (in a 
Bayesian setting). One high fidelity approach would be to abandon a correlation analysis 
in favor of a direct waveform template-based search for individual systems flM [15] .. 
along the lines of what has been proposed for detecting galactic binaries with a space 
based gravitational wave detector [IT]. The advantage of such an approach is that the 
signal model would accurately reflect the signals in the data, but the downside is that it 
greatly increases the size of the parameter space to be explored. We may find ourselves 
in a regime where each individual source lies below the detection threshold, while the 
combined signal may be detectable by some other less direct approach. A correlation 
analysis using a variant of (|8]), evaluated for several bright binaries with particular 
orientations and sky location, and with the frequency domain inner products restricted 
to sub-bands where one or two signals dominate, may be effective, but such an analysis 
introduces almost as many parameters as a multi-signal template based approach. One 
model that may find the sweet spot in the balance between fidelity and complexity 
introduces a single orientation parameter per bright source which helps to account for 
the anisotropy of the underlying signal. The orientation parameter is the angle between 
the actual signal direction and the signal direction used to construct the "correlation 
template" . To see how this is derived consider the filtered correlation function 

K,ji9, 0, 9t, 0t) = {rirj)ie, 0) A,-(^t, <Pt) , (20) 

with 

Pij{0T, 0t) = F+{eT, <pT)F+{eT, 0t) + {Ot, ct>T)F^ {Ot, 0t) . (21) 

The filter l3ij{6T, 0t) is the polarization averaged correlation function for a point source 
at sky location {9t, 4>t)- Note that sky average of this quantity is the standard HeUings- 
Downs correlation curve: 
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Averaging over Pulsar pairs separated by angle /i yields 

J Kij{9,(p,9T,(pT)S{cosfi-ni-nj)dQidQj = {H^)-f{fi,C) (23) 

where ( is the angle between the source direction [6, (p) and the filter direction {6t, (pr)- 
In the continuum limit, the standard p statistic is recovered by integrating the above 
expression over fj, and P = I (H^) lif^y C) d cos d cos (. The function 7(/x, () is plotted 
in Figure 5. Note that the matched filter 7(/i, 0) produces the largest correlation, and 
that using the sky averaged version of the filter {i.e. the average over () will degrade the 
sensitivity. In practice, since the source location is a priori unknown, it is not possible 
to parametrize the directional filter (21) by the angle (, and the search will have to 
be conducted using the parameters {9t,4>t)- But despite the parameterization being 
two-dimensional, the physical search is still one dimensional since a circle of points on 
the {9t, 4>t) sphere will yield exactly the same correlation curve 7(/i, (). 




Figure 5. The directional correlation function 7(/i, C)- 

For a black hole population the analysis could target the brightest black holes in 
each frequency band. For example, in the Bayesian formulation the correlation function 
to be used in the likelihood ( [l4| ) could be generalized to 

CM = E 0t) + SnXm, , (24) 

k 

where the S^{f) are localized to a particular frequency band. The optimal number of 
bands and their placement could be determined from the data using transdimensional 
Markov Chain Monte Carlo techniques. 
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6. Discussion 



We have shown that the standard cross-correlation analysis that was originally developed 
for isotropic, Gaussian backgrounds is capable of detecting the signals from individual 
black hole binaries, and by extension, the combined signal generated by an astrophysical 
population of binaries. We have also argued that the standard analysis will be sub- 
optimal in this case since the assumptions it makes about the signal are not valid, and 
we have suggested a number of approaches that may be more sensitive. We are currently 
exploring the relative performance of the various methods using simulated data from a 
variety of population synthesis models, and the results will be presented in a forthcoming 
pubhcation. 
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